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We report on the linear and nonlinear optical response of metamaterials evoked by first and 
second order multipoles. The analytical ground on which our approach bases permits for new 
insights into the functionality of metamaterials. For the sake of clarity we focus here on a key 
geometry, namely the split-ring resonator, although the introduced formalism can be applied to 
arbitrary structures. We derive the equations that describe linear and nonlinear light propagation 
where special emphasis is put on second harmonic generation. This contribution basically aims at 
stretching versatile and existing concepts to describe light propagation in nonlinear media towards 
the realm of metamaterials. 
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I. INTRODUCTION 

Metamaterials are attracting much scientific interest 
due to their unique optical properties in terms of their 
tailorable effective material response and the potential 
exciting applications [TJ [2l [3l H] . Most notably, the effec- 
tive magnetic response is an extraordinary property that 
distinguishes metamaterials from all naturally available 
media at optical frequencies. It has been shown that the 
effective material properties, such as the effective mag- 
netic and electric response, are determined by the res- 
onant excitation of plasmonic eigemodes, dependending 
on the specific metaatom geometry O E] . Frequently the 
structure and physics of these eigenmodes is investigated 
by calculating the electric fields or currents inside the 
metaatom by powerful numerical techniques. 
Here we present a more physical approach to charac- 
terize these eigenmodes by a multipole expansion up to 
the second order. This model permits for an analytical 
treatment of the light propagation in metamaterials at 
the expense that it constitutes merely a simplification 
of the real plasmon dynamics. Nevertheless, exactly for 
this reason the model provides unprecedented physical in- 
sights into the functionality of metamaterials which were 
inaccessible before by a pure numerical treatment. In 
the quasi-static approach for spherical or elliptical metal 
nanoparticles the fundamental plasmonic eigenmode is 
sufficiently described by an electric dipole. Metamate- 
rials, which in general exhibit more complex geometries 
with an involved carrier dynamics, require to account for 
higher order multipoles as e.g. the electric quadrupole 
and the magnetic dipole moment representing the natu- 
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ral extension towards a fully electrodynamic model [71 [8] . 
Furthermore the excitation of multipoles is frequently 
applied in order to explain the origin of the magnetic 
or the electric effective parameter response, phenomeno- 
logically. Here we present a coherent description corre- 
sponding to the way metamaterials are currently under- 
stood, namely by their multipole excitations. In par- 
ticular, the simultaneous consideration of the magnetic 
dipole and the electric quadrupole is required by nature, 
since both occur in the same order of the multipole ex- 
pansion. Though this has been extensively discussed in 
the literature [9l [lOl [Til HIH E] , the quadrupole mo- 
ment is frequently dropped. 

Besides the linear properties that can be covered by 
this expansion we present the extension of the multi- 
pole description towards the quadratic nonlinear opti- 
cal regime. Since the multipole expansion will be trun- 
cated beyond second order terms, we focus on the study 
of quadratic nonlinear effects associated with these sec- 
ond order multipoles. We mention that this proce- 
dure of introducing nonlinearity is known from the early 
works in nonlinear optics [15l [16] and is supported by 
several papers that observed multipole induced nonlin- 
ear optical effects in various plasmonic nanostructures 
[IZl[lll[Tl[20l|2ll|2a|2a[2lll^ 

Hence, the subject of this paper is twofold. At first we de- 
rive the required wave equations for the complete meta- 
material system and compare the results (the dispersion 
relation, and the effective material parameters) in the 
linear limit with numerical calculations for a representa- 
tive metaatom, i.e. the split-ring resonator (SRR). How- 
ever, the formalism we introduce is general and can be 
straightforwardly applied to other geometries. The SRR 
was only chosen because first experiments on the sec- 
ond harmonic generation were already reported in this 
structure; though in a configuration amenable for nano- 
fabrication where the induced magnetic dipole is non- 
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FIG. 1: (a) SRR metaatom and the intrinsic currents for the 
fundamental electric (black solid line) and magnetic (black 
dashed line) mode, (b) The associated auxiliary charge distri- 
bution (red points) with predefined degrees of freedom (black 



radiating [26} [27] . To become specific in the second part 
second harmonic generation (SHG) in a SRR is studied 
in detail. In order to predict the second harmonic gen- 
eration by purely analytical means, we finally apply the 
undepleted pump approximation (UDPA) and derive the 
associated equations within the multipole model. We 
observe an excellent agreement with the numerically de- 
rived solutions which justifies the application of this ap- 
proximation for further predictions. 
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where the equations of motion r/(r,t) for each charge qi 
contain all information about the plasmonic eigenmodes. 
The parameter r] accounts for the dipole density. Even 
at this early stage it can be seen, that if r/(r, t) (x E(r, t) 
second order multipoles will immediately evoke nonlinear 
contributions, since they involve terms (x r^(r, t). For the 
carrier configuration proposed here the associated terms 
ri{r,t) are 



(-xo,^o,0), rj" = (-xo -^1,^0 -^1,0), (3) 
(xo, yo, 0), = {xq - 6, ^o, 0), 
(-xo, -^0, 0), r3 = {-xo - 6, -yo + 6,0), 
(xo, -yo, 0), = (xq - 6, -^o, 0). 



II. THE NONLINEAR WAVE EQUATION 



Our investigation starts with the wave equation incor- 
porating multipoles up to second order [13] 



AE(r,t) = 



(1) 



^2 ^ 

^oV • ^Q(r, t) + /ioV X ^M(r, t). 



Now, in terms of the plasmonic eigenmodes of inter- 
est the respective metaatom has to be mapped onto 
the point multipoles: electric dipoles P(r,t), magnetic 
dipoles M(r, t) and electric quadrupoles Q(r, t). In order 
to observe both, an electric and a magnetic response, the 
SRR is uprightly oriented [6 , see FIG.Q. To cover the 
fundamental electric and magnetic modes [5 , sketched 
by the black solid and dashed lines in FIG.([l^), respec- 
tively, four auxiliary positive and negative charges with 
predefined spatial degrees of freedom are required as in- 
dicated in FIG.([l]3). With the knowledge about these 
directional constraints for the carrier dynamics the mi- 
croscopic definitions of the multipole moments can be 
applied as 



In Eq.(|3| the superscripts ± denote whether the po- 
sition vector is associated with a positive or a negative 
charge, while the argument (r, t) is suppressed for all r^"^. 
Furthermore, all positive carriers are fixed at the posi- 
tions ±xo, ±^0, while negative carriers are allowed to os- 
cillate around these sites, described by = 6,2(1*, t). 
The distinction between the carrier oscillations in both 
SRR wires is vital for realizing the two plasmonic eigen- 
modes [FIG.([l^)]. These carrier oscillations evoked by an 
external electromagnetic field and intrinsic Coulomb in- 
teractions may be described by a set of coupled oscillator 
equations as 



m 
m 



E^{y^yo,t), (4) 

Ex{y - yo^t)' 



In Eqs.Q 7 represents the damping, loq the eigenfre- 
quency while a describes the coupling strength between 
the carriers in the SRR arms. The physical origin of 
this coupling is the Coulomb interaction of carriers in 
the horizontal SRR arms excited by an electric field par- 
allel to the arms and the carriers in the vertical arm that 
are excited by the local fields of the horizontally oscillat- 
ing charges, producing a current inside the entire SRR. 
Moreover this coupling between two identical oscillators 
results in a splitting into symmetric and anti-symmetric 
oscillation modes. Substituting Eqs.Q into Eqs.([2| the 
remaining multipole moments are obtained as 
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P(r, t) = e,2r?q(ei + 6) + 63,^(^1 " 6), (5) 
M{r,t) = -e,^(xo + 2yo)^(Ci-6), 
Q.y{r,t) = |?(a-C2)[(2yo-a:o)-(Ci+6)]- 

From Eqs.([5| it can be deduced that all multipoles de- 
pend either on the sum or the difference of and ^2- 
Especially for the symmetric carrier oscillation in the 
SRR arms (^i = ^2) all second-order moments vanish 
and only two identical electric dipoles parallel to the SRR 
arms remain (symmetric mode). In turn, an antisym- 
metric oscillation (^1 = —^2) excites both second order 
multipoles and a longitudinal electric dipole only since 
the electric dipoles in the top and bottom SRR arms 
are canceling each other (antisymmetric mode). Thus, 
the charge alignment chosen meets all requirements to 
describe the desired plasmonic eigenmodes and their dy- 
namics is determined. Decomposing the electric field into 
plane waves at the fundamental (FF) and the second har- 
monic frequency (SH) as E(r, t) = exEx{y^t) with 



the solutions to the oscillator equations [Eqs.Q] read as 
a ± 6 = ^±i?„e^('^(")^-'^*) (7) 

where the amplitudes are given by 

C = 2xicos(/c(c^)^o), (8) 
C = 2ix;;sin(A:(cj)?/o), 

with the introduced quasi-susceptibility 
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m uJq — uj"^ — zjuj ± a 

The respective equations for the SH field follow by substi- 
tuting u by 2uj in Eqs.([8|. In contrast to ordinary electric 
dipole interaction we observe a frequency splitting in the 
quasi-susceptibility evoked by the two frequency de- 
generated eigenmodes. Now, Eqs.([5]-[9| can be inserted 
into Eq. ([T]) yielding a nonlinear eigenvalue equation with 
second order nonlinear source terms 
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where the following abbreviations have been used for the 
linear 
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and the nonlinear multipole source terms 
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The exact solution to this eigenvalue equation 
would result in a nonlinear dispersion relation with 
k{2u, E^^2uj), k{uj, E^^2uj) and a fixed ratio E^/E2uj' The 
left-hand side of Eqs.(ITo| contains a part well-known 
from dipole interaction Xplo) but in addition contributions 
which stem from the second-order multipole response 
{(luji^uj)' Additionally the quadrupole moment causes 
a nonlinear term on the right-hand side. Interestingly, 
in this model the nonlinear response of the magnetic 



dipole produces no nonlinear contributions which is 
supported by rigorous simulations for a corresponding 
SRR configuration There the magnetic nonlinear 

contributions have been shown to be much smaller in 
comparison to a convective electric current [29 which 
is equivalent to the quadrupole contribution in our 
approach [37]. Furthermore, it is mentioned that in 
contrast to usual second order nonlinear optics [30 here 
the nonlinear source term incorporates the first spatial 
derivative, induced by the quadrupole moment. 



III. LINEAR OPTICAL PROPERTIES - 
EFFECTIVE MATERIAL PARAMETERS 

In order to validate the predictions of the model we 
start with the investigation of the linear properties. To 
this end the nonlinear source terms [Eqs.(12)] have been 
dropped which yields two decoupled linear eigenvalue 
equations for the FF and the SH wave. The obtained 
linear wave equations describe the field propagation in 
an effective medium determined by its multipolar contri- 
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FIG. 2: (a) Comparison of the dispersion relation - numerical 
model (dashed), analytical model (solid); (b) Single meta- 
material layer (gold SRR) used for the numerical simulation 
(dimensions in nm); spectral dependence of effective permit- 
tivity (c) and effective permeability (d) from the multipole 
model; spectral dependence of effective permittivity (e) and 
effective permeability (f) from the numerical simulations. 



butions. Thus the correlated wave vector k{uj)^ i.e. the 
dispersion relation, represents a self consistent solution 
and contains all physical information to describe light 
propagation in the respective metamaterial. This can be 
understood in complete analogy to the exact parameter 
retrieval which assumes also a homogeneous wave propa- 
gation inside a metamaterial slab. By neglecting the the 
nonlinear source terms qu;2u,-uj and q2ou;u,uj in Eq.(lO) 
we get the linear wave equation 



d 



xe 



dy_ 



For plane waves their solution yields the linear dispersion 
relation 



ik{uj) 



CqUJ 



. (14) 



Together with the multipole definitions [Eqs.(pT])] this 
implicit equation becomes 



k\u) = -^[1 



, ^VQxi cos{k{uj)yo) 
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(15) 



which has to be solved numerically. However in the sub- 
wavelength limit yo < X the trigonometric functions can 



be expanded 
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Thus the linear dispersion relation and hence the effective 
refractive index may be explicitly written down as 
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The remaining unknown parameters of the introduced 
formalism can be obtained by fitting this dispersion re- 
lation to the numerically calculated one for the present 
gold SRR geometry, shown in FIG.([2|l) [3T1. As a result 
the two first resonances [denoted by (1) and (2)], which 
are the fundamental magnetic (1) and the fundamental 
electric (2) mode, are nicely covered by our dispersion 
relation. It can be clearly seen that the next higher mag- 
netic mode [resonance (3)] is not covered by the multi- 
pole description. This would require to consider a more 
complex carrier dynamic as it can be deduced from nu- 
merical simulations [5 . Moreover, the effective param- 
eters (ceff, /ieff) follow without any further fitting [13] 
and show a good agreement in comparison with the cor- 
responding numerically determined values FIG.([2]3-f). In 
order to determine the effectice permittivity and perme- 
ability too we have to consider the constitutive relations 
[32 , 33 , by which we introduce both quantities. We start 
with the electric permittivity 
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In complete analogy we derive the effective magnetic per- 
meability and begin with 



1 



M, 



H(r,t) = ^B(r,t) 
Mo 

H,{y,t) = ^B,{y,t)-M,{y,t), 
/^o 

— — m^E^. 
Mo 



(19) 



At this point it becomes evident that the effective mag- 
netic response is merely evoked by the electric field [34 . 
This is a consequence of the underlying model, because 
the complete carrier dynamics is driven by the electric 
field only. Within the multipole expansion this carrier 
dynamics, i.e. the plasmonic eigemodes are transformed 
into multipole moments including the magnetic dipole 
moment. Thus, all moments, even the magnetic moment. 
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are a mere consequence of the electric field interaction 
with microscopic carriers. We mention that this is oc- 
casionally erroneously interpreted, by an argumentation 
that the magnetic moment is induced by the magnetic 
field. Our approach shows that all effects, also magnetic 
ones, may be sufficiently well described by an electric 
field interaction only. Nevertheless, to define the mag- 
netic permeability it is necessary to translate the electric 
field into a magnetic field which yields 



1 

Mo 



1 + m,. 



1 



1 + m, 



(20) 



In contrast to the retrieval of the systems parameters we 
have fitted the disperison relation to the model, but in 
turn this can be performed equivalently for the effective 
material parameters eeff and /ieff, yielding the same phys- 
ical information with respect to the required quantities. 



IV. NONLINEAR OPTICAL PROPERTIES - 
SECOND HARMONIC GENERATION 

To study the nonlinear behavior induced by the fun- 
damental modes of the present metaatom we resort to 



the linear dispersion relation and treat the nonlinearity 
as perturbation rather than solving Eqs.(lO) exactly. As 
usual we rely on the slowly varying envelope approxima- 
tion (SVEA) [30^. Within this approximation the fast 
spatial oscillation ex.-p[ik{uj)y] is separated from a slowly 
varying amplitude A{y), which contains all information 
about the generation and depletion of the fundamen- 
tal {E^ = A^{y) ex.-p[ik{u;)y]) and the second harmonic 
{E2u; = A2uj{y)exp[ik{2ij)y]). 



A. Exact numerical solution 

At first the solution to the wave equations incorporat- 
ing the SVEA ansatz has been performed numerically. 
Therefore we simplify this system by introducing the fol- 
lowing substitutions 



Mo (' 
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(21) 



Now the eigenvalue equations take the following form 
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(22) 



The solution to these equations are two coupled non- 
linear dispersion relations, one for the fundamental and 
one for the second harmonic wave, which depend on both 
fields. In order to avoid the solution of this involved sys- 
tem we treat the nonlinearity as a perturbation and re- 
sort to the linear dispersion relation. To study the effect 



of the nonlinear source terms we apply the slowly vary- 
ing envelope approximation where the linear fields are 
weighted by a slowly varying amplitude functions A(y). 
Replacing the constant amplitudes E^^2uj by A^^2uj{y) 
we obtain upon substitution and upon neglection of the 
second-order derivatives (SVEA) 



^AUy) = -^||ff5^ (^A2Uy)^AUyy^AUyy^A2Uy)^^{k{2u)-k^^ 



dy 

^ ^i{k{2uj)-k{u;)^ -k{u))y ^ 



I 

This final system has been solved numerically [FIG.(|3])]. The FF wave evolution E^{y) is determined 
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FIG. 3: Evolution of normalized electric field intensity for 
the FF (a) and the SH (b) as a function of the wavenumber. 
The red lines indicate the real (dashed) and the imaginary 
part (solid) of the linear dispersion relation, (c), (d) The cor- 
responding results for the un-depleted pump approximation 
(UDPA). 



by the two eigenmode resonances (both indicated by the 
resonances in the red lined dispersion relation), where a 
strong damping is observed. For frequencies out of the 
spectral domain of these resonances the FF wave propa- 
gates without excessive losses, as expected. For the SH 
wave a strong contribution at the fundamental magnetic 
and electric resonance can be observed. In these calcu- 
lations the SHG signal originating from the electric reso- 
nance around u = 6400 cm~^ seems to be much stronger 
than in the spectral vicinity of the magnetic resonance 
(p = 4000 cm~^). This originates from the strong damp- 
ing of the SHG wave for the magnetic resonance (which 
propagates ai u ^ 8000 cm~^), because in this spectral 
domain an enhanced damping occurs due to the presence 
of the electric resonance. This changes dramatically for 
the SH wave induced by the electric mode, since at the 
SH frequency the imaginary part of k{2uj) is close to zero. 
Thus the second harmonic originating from the electric 
resonance propagates almost without damping. 



B. Un-depleted pump approximation 

In order to double-check the results and to take the 
weak conversion efficiency into account the undepleted 
pump approximation (UDPA) has been applied [30] . 
Within this approximation the fundamental wave (the 
pump) remains unaffected by the generated second har- 
monic wave. This can be expressed by setting tlJou;2ou,-u in 



the first equation as well as the first derivative of A^{y) 
in the second equation to zero. This results in 



d_ 
dy 
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1 



^i{2k{uj)-k{2uj) 



(24) 



Comparing the numerically determined electric field for 
the fundamental and the second harmonic wave to those 
of the undepleted pump approximation [see FIG.([3|] one 
can clearly deduce that the UDPA describes the propaga- 
tion for both waves almost exactly. Furthermore Eq.(24) 



permits to calculate analytically the conversion efficiency 
from fundamental to second harmonic intensity, e.g. for 
a slab consisting of a single layer of SRR's. It is impor- 
tant to note that for this purpose only parameters are 
required, that are fixed by comparison with the linear 
effective material response. In passing we comment that 
such a procedure; the determination of the nonlinear ma- 
terial properties based on the linear material parameters, 
is known as Miller delta a well established rule in non- 
linear optics [35l |36J . 



V. SUMMARY 

In summary, we have presented a self-consistent physi- 
cal model that permits to describe the linear response of 
metaatom geometries by their intrinsic plasmonic eigen- 
modes. The occurring specific carrier dynamics have 
been mimicked by an auxiliary carrier alignment interact- 
ing with the incident radiation. The knowledge of these 
charge oscillations allows the application of the multi- 
pole expansion which provides the eigenvalue equation 
for electromagnetic waves propagating in such a com- 
posite metamaterial. Moreover we have shown that the 
specific convective carrier oscillations together with the 
quadrupole moment inherently introduce nonlinear ma- 
terial interactions. Considering the SHG process, our 
calculations show the expected enhanced signal both for 
the electric and the magnetic resonance and provide a 
microscopical and physical understanding of them. For 
further investigations it is important that the nonlinear 
response can be determined only from knowing the linear 
response, which is accessible by comparing the dispersion 
relation or any other effective material property to the in- 
troduced multipole model. 
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